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O : ABSTRACT 

o . 

(N ■ 

• We describe a new method for analyzing gravitational lens images, for the case 
pL| ' where the source light distribution is pixelized. The method is suitable for high 
OO ' resolution, high 5"/^^ data of a multiply-imaged extended source. For a given 
^ ■ mass distribution, we show that the step of inverting the image to obtain the 
^ ■ deconvolved pixelized source light distribution, and the uncertainties, is a hnear 

' one. This means that the only parameters of the non-linear problem are those 

^ ■ required to model the mass distribution. This greatly simplifies the search for a 

^ ■ min.— fit to the data and speeds up the inversion. The method is extended in 

^ ' a straightforward way to include linear regularization. We apply the method to 

O ■ simulated Einstein ring images and demonstrate the effectiveness of the inversion 

• for both the unregularized and regularized cases. 

2 ' Subject headings: gravitational lensing 

c/3 '. 



1. Introduction 

This paper is concerned with the problem of inversion of a gravitationally-lensed image 
of an extended source, i.e. a galaxy rather than a star or quasar. This problem is interesting 
because lensed images of extended sources provide more information than images of point 
sources, and so potentially are more useful for determining the mass profiles in galaxies and 
clusters of galaxies. Also, because of the magnification, one can measure structure in the 
light profile of the source at enhanced resolution. In this paper we show how this problem can 
be separated in a natural way into linear and non-linear dimensions, and how this provides 
a number of advantages. 

In this introduction we review solutions to the inversion problem and introduce some of 
the terminology used in the remainder of the paper. In all the solutions described here the 
mass in the lens is parameterized. Nevertheless the analysis applies equally to a pixelized 
mass distribution. 
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When presented with the lensed image of an extended source, the unknowns to solve for 
are the source hght profile and the lens mass profile, and the uncertainties in these quantities. 
One approach to this problem, suggested by Kayser and Schramm (1988), uses the fact that 
regions of the source that arc multiply-imaged have the same surface-brightness. For a 
trial mass distribution, the method traces image pixels to the source plane where the counts 
in different image pixels mapping to the same source pixel are compared. The solution 
for the mass is obtained by minimizing the dispersion in the image pixel counts for such 
multiply-imaged source pixels. Kochanek et al. (1989) successfully applied this approach 
to the inversion of the radio Einstein ring MGl 131+0456. The algorithm was refined by 
WaUington, Kochanek, and Koo (1995) who apphed it to the triply-imaged giant arc in the 
galaxy cluster CI 0024+1654. 

The main shortcoming of this approach is that it does not deal with the image point- 
spread-function (psf). If psf smearing of the image (either instrumental or atmospheric) 
is significant, the light profile of the source is not correctly recovered by backward tracing 
the image, even if the mass distribution is exactly known. To deal with the psf, a forward 
approach is needed i.e. one chooses a model for the source light profile (parameterized or 
pixelized), and a model for the mass (parameterized or pixelized), forms the image, convolves 
it with the psf, and compares it to the actual image, adjusting the source and lens models 
to minimize a merit function e.g. x^- 

An argument for choosing to parameterize rather than pixclizc the source light profile 
is that it forces the solution to be smooth. Nevertheless, the source light profile may be 
complex. This is true, for example, in the cases of MG1131+0456 and CI 0024+1654, cited 
above. A large number of parameters might be required to provide a satisfactory description. 
Without clues to the character of the source it is extremely difficult to select the best 
parameterization i.e. the one which provides a satisfactory fit with the smallest number 
of parameters. In the most extreme example Tyson, Kochanski, & dell' Antonio (1998) used 
232 parameters to model the source light distribution of the galaxy lensed by the cluster 
C10024+1654. 

If the source light profile is complex it is natural to consider pixelizing the source, 
i.e. the counts in each pixel is a free parameter. This removes the difficulty in finding a 
good parameterization for the source, and thereby avoids any bias in the fitted mass profile 
resulting from a poor choice. On the other hand, due to the deconvolution, and because 
the pixels are independent, the solution can be noisy. It is possible to achieve a smooth 
pixelized solution by adding a suitable 'regularizing' term to the merit function. This term, if 
minimized on its own, would produce a smooth source light profile. By adding this term to 
the final solution involves a balance between obtaining the best fit to the image (minimizing 
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X^), and obtaining a smooth source solution (minimizing the regularizing term). Wallington, 
Kochanck, and Narayan (1996) apply this approach to the case of the radio Einstein ring 
MG 1654+134. They use a maximum entropy approach i.e. the regularizing term to be 
minimized is the negative of the entropy, the negentropy. Labeling the counts in source pixel 
i by Si, the source plane negentropy is ^^Siln(si), and the merit function they minimize is 



(here we have followed the notation in Press et al., 2001, §18.7). The purpose of the multiplier 
A is to give more or less weight to the negentropy term. 

The inversion proceeds as follows: For a fixed value of A, the solution is determined 
by searching through the multi-parameter space for the minimum of the merit function. 
The number of dimensions of the parameter space to search is the sum of the number of 
source pixels and the number of parameters used for the mass. This search is most efficiently 
achieved with a pair of nested cycles. The inner (source) cycle searches for the best source 
light profile for a fixed mass profile. The outer (mass) cycle adjusts the mass profile. Outside 
this cycle is a third (A) cycle where the multiplier is adjusted. Because the negentropy term 
acts to smooth the source, as A increases, the xfm term also increases, i.e. the fit becomes 
worse. The principle for reaching the final solution (e.g. Press et al., 2001, §18.4) is to start 
with A large, then progressively reduce the weight of the regularizing term until the xlm 
becomes satisfactory. In other words the solution has the smoothest source that provides a 
satisfactory fit to the image. 'Satisfactory' is usually interpreted as reaching the criterion 
for the for the image xtm — "^^'^(xLi) + '^(xfm)- With three nested cycles, the A, mass, 
and source cycles, the routine can be slow. 

In this paper we describe a new technique which we suggest simplifies and clarifies 
the problem in a number of ways. In purely formal terms, the method is very similar to 
the maximum entropy method of Wallington et al.: Algebraically, we simply replace the 
negentropy term in the merit function (1) with a linear regularization term. However, the 
insight we bring is to show that for a fixed mass distribution, the minimization of the merit 
function is now a linear problem i.e. can be solved by matrix inversion. In other words 
the source cycle - the innermost of the three minimization cycles - is eliminated. This 
has major benefits. In the first place the inversion is much quicker, thereby allowing a more 
thorough search for the best fit mass model. At the same time, the uncertainty of identifying 
the true minimum has been removed. The method also greatly simplifies calculation of the 
uncertainties, as we show below. More generally, the formalism clarifies the essence of the 
problem: The source parameters are linear dimensions and the mass parameters are non- 
linear dimensions. For this reason we call the method 'semi-linear'. 




(1) 



-4- 



At this point it is worth noting that, because of magnification and multiple imaging, the 
number of constraints to the solution can be much greater than the number of parameters 
to solve for. In this respect the lens inversion problem differs from many inversion problems 
encountered in astronomy (for example image deconvolution). We find, as a consequence, 
that in many circumstances the regularization term can be removed altogether. So the 
A cycle is also eliminated. The merit function is then just xf^, and this is our starting 
point for the presentation of the theory. For a fixed mass profile, the pixehzed source light 
distribution that produces the min.— x?^ fit is obtained by linear inversion. The mass profile 
is then adjusted to find the minimum of these min.— fits. The advantage, besides speed 
(only the mass cycle remains), is that the solution is unbiased, since there is no smoothing 
of the source. 

The outline of the remainder of the paper is as follows. In §2 we explain the basic theory, 
demonstrating that for a fixed mass profile, the problem of obtaining the source light profile 
by x^ minimization in the image plane is a linear one, and obtaining the covariance matrix 
for the counts in the source pixels. We then extend the basic theory to include a linear regu- 
larization term. In §3 we apply the method to a realistic problem, assessing the performance 
for different psf widths, and different source pixel sizes, with and without linear regulariza- 
tion. In §4 we provide a summary of the main points, together with recommendations for 
applying the method. 

2. Theory 

In this section we present the theory of semi-linear inversion, firstly without regular- 
ization (§2.1), and then with regularization (§2.2). In each sub-section we begin with the 
case where the mass is fixed, and then treat the general case, minimizing also on the mass 
parameters. 

2.1. Semi— lineeir inversion without regularization 

2.1.1. Fixed mass: Eliminating the source cycle 

Without any regularizing term, the merit function is simply G = Xim- The basic problem 
is to find the counts in the source pixels that, for a given mass distribution, minimize the 
merit function G, i.e. give the best fit to the observed image. Pixels in the source plane are 
labeled i = 1,1. There is no restriction on how the source plane is tessellated. In principle, 
the pixels could change in both size and shape across the source region, which itself could be 
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of any shape. Pixels in the image plane are labeled j — 1, J. It is assumed in the following 
that the image pixels include counts from the image of the lensed source only i.e. the images 
of any Icnsing galaxies, and the mean sky count, have been subtracted. Also we suppose 
that the data in each image pixel are independent i.e. are characterized by the counts dj, 
and dispersion aj, with no covariance between pixels (appropriate for CCD data). 

The inversion proceeds as follows: Choose the mass model parameters, then, for each 
source pixel i, in turn, form the image for unit counts (surface brightness) by appropriate 
ray tracing and convolution with the known point spread function i.e. compute the counts 
in the ith image fij,j = 1, J. The problem now is to combine these / images with scalings 
Si,i = 1,1, to minimize G. These scalings are the deconvolved intrinsic source surface- 
brightness distribution. 

The problem is of a standard type. The merit function is written 

2 

(2) 
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Aim 



E 
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Minimizing with respect to each of the source terms we have a set of I simultaneous 
equations of the form 



fij X]fe=l ■'^kfkj fijdj 



-I 



(3) 



symmetric 7x7 matrix, with elements Fjfe = X]/=i fijfkj/'^j- Finally D is a column matrix 
of length 7 containing the elements Dj = '^i=i fijdj /cr]- 



where the reason for the factor | will soon become clear. These equations may be written 
in matrix form 

FS = D. (4) 

Here S is a column matrix of length 7 containing the elements Sj, to be solved for. F is a 

.j=i fijdj /(^l 

The solution for the counts in the source pixels is then simply obtained by matrix 
inversion 

S = F-^D (5) 

thus eliminating the source cycle. 

The solution for the errors has a particularly simple form. We seek the covariance matrix 
for the source pixels. Noting that 
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we see that the matrix F is one half times the Hessian matrix of x?^, which is to say that F 
is the curvature matrix of the problem (Press et al., 2001, §§15.4, 15.5) - this was the reason 
for using the factor | in equation (3). We now show that the matrix C = F~^ is the required 
covariance matrix of S. 

For independent image pixels, the covariance between source pixels i and k is given by 
Using equation (5) this becomes 

J / / ■ ^ / • 

'^ik ^J2^^Y1 ;i ^ • 

j=l 1=1 j m=l j 

Multiplying this out gives 

4 = Cik (9) 

as required. 

We see that for the case of fixed mass, the covariance matrix for the source pixel counts 
is computed in the inversion step, without the need for further calculation. We shall refer 
to this I X I matrix as the source covariance matrix hereafter. Even though it is not the 
complete solution for the source pixel errors (because the mass parameters have been fixed), 
the source covariance matrix is extremely useful, for example, in exploring different mass 
models and pixelizations (§3). 

It is worth noting that the scmi-lincar inversion solution, cither with or without reg- 
ularization, differs in character from the maximum entropy solution. With the semi-linear 
method the counts in any source pixel are unbounded, so the best-fit value could be neg- 
ative, since some image pixels contain negative counts (i.e. are below mean sky). With 
the maximum entropy method all source counts must be positive. The semi-linear method 
provides the best estimate of the counts in a source pixel, and the solution is satisfactory 
provided the result is consistent with being positive. If the counts in any source pixel are 
significantly negative (e.g. < — 4o") this would indicate a bad mass model. This possibility 
of testing the suitability of the mass model with a source-plane statistic can be viewed as 
an extra positive feature of the semi-linear method. 
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2.1.2. Mass cycle 

The full solution proceeds by searching through the mass-distribution parameter space, 
at each point minimizing xlm by linear inversion, to find the smallest of these min.— 
values, the global minimum. Because the number of dimensions of the parameter space for 
the non-linear search has been greatly reduced, it is now a much simpler problem to locate 
the true minimum securely. 

The solution for the errors is more complicated than above, since we have added in the 
non-hnear mass dimensions. If there are L parameters that describe the mass, labeled m;, 
we need to form the (/ + L) x {I + L) (symmetric) curvature matrix of the problem. But note 
that the majority of the terms, the 7x7 terms |gf^-, have already been computed and are 
the elements of the matrix F at the global minimum. The remaining terms, the L rows (and 
columns) of terms such as ^ „ • and hjr-^, can be filled in by simple measurement 
of the shape of the x^„j surface about the global minimum. The covariance matrix for the 
mass and source parameters is the inverse of this curvature matrix. We shall refer to this 
(7 + L) X (7 + L) matrix as the full covariance matrix hereafter. 



2.2. Semi— linear inversion with regular izat ion 

2.2.1. Fixed mass: Eliminating the source cycle 

The possibility of replacing the negentropy term in equation (1) by a term (a linear 
regularization term) which preserves the linearity of the min.— approach is made evident 
by the linearity of equation (3) with respect to the source parameters. Clearly we can form 
a merit function by adding to any term Gl which is a linear combination of terms siSk 

GL = ^aikSiSk (10) 

i,k 

since the partial differentials of these additional terms will also be linear. One example of a 
linear regularization term is Gl = Yl!i=i ^l- The choice of Gl is discussed below. 

Writing the merit function generally as 

G = xL + AGi (11) 

then, following through the same analysis as in §2.1.1, the solution for the counts in the 
source pixels can be written 

S = [F + AH]-^D. (12) 
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We call the matrix H the regularization matrix. The elements of H are 




matrix. 

The form of Gl should be chosen to penalize noisy solutions. The choice Gl = Y^\=i sf, 



termed "zeroth-order" regularization in the literature, is one attempt to achieve this. Other 
widely-used linear regularization terms include gradient and curvature forms. These three 
regularization forms correspond, loosely speaking, to the prejudice that the source light 
profile is, respectively, approximately zero, constant, or planar (see Press et al., 2001, §18.5, 
for a detailed account of the theory of linear regularization and its implementation). In 
practice, if A is not too large, all three terms serve to smooth the source in a rather similar 
way, and there is little to distinguish between the solutions. 

The gradient and curvature forms consider the differences between counts in neighboring 
pixels. Until now we have used a one-dimensional numbering scheme for the source pixels. 
In this case, since we need to take account of the relative spatial locations of pixels in the 
source plane we use coordinates x,y. The simplest gradient term is 



Another form uses [sx,y — 0-5{sx+i,y + Sx,y+i)]'^- In forming the sum it is necessary to translate 
the indices x,y to the index i, and then equation (13) is used to form the matrix H. Note 
that zeroth-order regularization is computationally by far the simplest method, since it does 
not involve this step of translation of indices. 

The simplest curvature form is 



Another form uses [sx,y - 0.25{sx-i,y + Sx+i,y + Sx,y-i + Sx,y+i)]'^. 

The source covariance matrix for the regularized case is fortunately only slightly more 
difficult to compute than for the unregularized case. Writing R = [F + AH]~^, and following 
the same line of reasoning as in §2.1.1, we obtain the analogous equation to equation (8) 
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Multiplying out we obtain 



I 




kl- 



(17) 



1=1 



2.2.2. 



Mass cycle 



The procedure for the full solution is the same as for the unregularized case. One 
searches through the mass-distribution parameter space, at each point minimizing G by 
linear inversion, to find the smallest of these min.—G values, the global minimum. In the 
regularized case there is no simple solution for the full covariance matrix however. In the 
unregularized case, we were able to use the fact that the inverse of the full curvature matrix 
is the full covariance matrix. But in the regularized case this is no longer true since we are 
minimizing G = xlm + ^l- Instead, an alternative approach must be used, for example, a 
Monte Carlo method which inverts an ensemble of realisations of the image by adding noise 
to the original image. 



In this section we apply the semi-linear inversion method to a realistic test problem. 
To validate the linear inversion step, we begin with the case of fixed mass. We quantify the 
effectiveness of the method under variations of the image S/N, psf width, and source pixel 
size, for both the unregularized and regularized cases. We then present an analysis of the 
full problem, allowing the mass parameters to vary. Finally we debate the advantages of the 
unregularized and regularized approaches, for different practical applications. 



To make the computations more useful we have based our investigation on a realistic 
simulation of a deep image of an Einstein ring gravitational lens system observed with the 
Advanced Camera for Surveys (ACS) aboard HST. The camera has a pixel size of 0.05". We 
have used cosmological parameters = 0.3, Qa = 0.7. The lens is placed at z — 0.3 and 
the source lies at 2; = 3.0. 

Figure 1 illustrates the test problem. The lens (not shown) is modelled as a singular 
isothermal ellipsoid with one-dimensional velocity dispersion 260 km s~^ and ellipticity e = 



3. 



Simulations 



3.1. Test problem 
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1 — b/a — 0.25. The semi-major axis of the lens is ahgned at 40° counterclockwise from 
the vertical. The Einstein angle is 9e — 4:T^cr^Dds/{c^Ds) = 1-58". The source, shown in 
the top left panel, is contained within a square of size 0.75" and is modelled as two circular 
sources of Gaussian profile, binned in 0.05" pixels, the same as the image pixel size. The 
peak surface brightness of each source is 1.0, in arbitrary units. One source lies inside the 
inner caustic, while the second source straddles the inner caustic. This source configuration 
resembles that inferred for the gravitational lens 0047—2808 (Wayth et al., 2003). To create a 
realistic ACS simulation the image was formed by ray tracing, then convolved with the point 
spread function, and noise added (in the manner described in the following paragraph). For 
simplicity we modelled the psf as a Gaussian, and chose FWHM 0.08" which is the resolution 
of a diffraction-limited telescope of the same diameter as HST, at a wavelength A = 800nm. 
Because of the slight undersampling, the convolution is made on a sub-pixelized grid and 
then binned up to the full pixel size. 

The data pixels used for the inversion were the 3626 pixels within the annulus in the 
image plane marked in the figure. This annulus is defined by the region covered by imaging 
the entire source plane. ^ An important point to note is that the analysis region must at least 
cover this annulus, otherwise the counts in some source pixels will be unconstrained and the 
inversion will fail. A larger region may be used, but if it becomes too large the usefulness of 
the statistic is diminished, as then most of the pixels are in the background. The final 
step in the simulation is to apply uniform Gaussian random noise over the image plane. The 
noise level is defined in terms of the total S/N^^ integrated over the annulus. The same 
noise realisation was used for all the simulations, but scaled in order to vary S/Ni^. 

The upper middle panel shows the final simulated ACS image. This image, with source 
pixel size 0.05", S/Ni^ — 60, and psf FWHM= 0.08", is the reference test problem to invert. 
We later vary these three parameters. The parameters of the different models we have run 
are listed in Table 1. Col. (1) provides the simulation number. The reference problem is 
numbered 1. Col. (2) gives the source pixel size in arcsec, col. (3) the summed S/N in the 
image, and col. (4) the psf FWHM in arcsec. Col. (5) is a label U or R depending on whether 
the inversion was unregularized or regularized. Then, in the case of regularized inversion, 
col. (6) provides the degree of regularization, quantified by the parameter iV^, which is the 
increase of Ximi^) from the minimum in units of the standard deviation o"(x^^(z/)). Recalling 
the discussion in §1, a value Nx — 1 in this column corresponds to the usual criterion for 
the maximum allowed degree of regularization. The other columns are explained in the next 
section. 



The region of the central image should also be included for non-singular mass models. 
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Fig. 1. — This plot shows the unregularized solution for the reference problem, line (1), Table 
1. The source plane, top left panel and bottom row, is 0.75" x 0.75" with 0.05" pixels, and is 
centered on the optic axis. The source comprises two circular Gaussian components and is 
shown top left. Also marked is the line of the inner caustic for the isothermal ellipsoid lens. 
The image, convolved with the psf, FWHM 0.08", and with noise added, is shown upper 
middle. The image pixel size is 0.05" and the image box size is 5.0" x 5.0". The lower left 
panel is the source light distribution reconstructed from the image by semi-linear inversion 
without regularization. The upper right panel is the image of this source, convolved with the 
psf. The lower right panel displays the la uncertainty for the source pixels, and the lower 
middle panel is the source S/N image. The dotted square is the region over which S/Ngo 
is measured. In this and the following two figures counts in pixels in both the image and 
source plane are in units of surface brightness. 
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(1) 

sim. 
no. 


(2) 
source 

pix. size 
// 


(3) 

S/Nim 


(4) 
psf 

FWHM 
// 


(5) 
U/R 


(6) 
Nx 


(7) 
X^M 


(8) 

Xloi^) 


(9) 


(10) 
|As/c7| 
max;. 


(11) 


(12) 
note 


1 


0.050 


60.0 


0.08 


U 




0.956 ±0.024 


1.088 ± 0.094 


79.9 


2.80 


0.082 


Fig. 1 


2 


0.050 


30.0 


0.08 


u 




0.956 ± 0.024 


1.088 ± 0.094 


40.6 


2.80 


0.163 




3 


0.050 


60.0 


0.00 


u 




0.958 ±0.024 


1.052 ±0.094 


111.2 


2.64 


0.040 




4 


0.050 


60.0 


0.16 


u 




0.956 ±0.024 


1.090 ± 0.094 


24.6 


2.73 


0.252 




5 


0.025 


60.0 


0.08 


u 




0.952 ±0.027 


1.003 ±0.047 


20.2 


2.84 


0.634 


Fig. 2 


6 
7 
8 
9 
10 


0.050 
0.050 
0.025 
0.025 
0.025 


60.0 
60.0 
60.0 
60.0 
60.0 


0.08 
0.08 
0.08 
0.08 
0.08 


R 
R 
R 
R 
R 


1 

2 
1 
3 
5 


0.980 ± 0.024 
1.004 ±0.024 
0.979 ±0.027 
1.033 ±0.027 
1.088 ±0.027 


1.138 ±0.094 
1.461 ± 0.094 
1.003 ±0.047 

1.003 ±0.047 

1.004 ±0.047 


111.8 
120.8 
64.5 
86.6 
98.2 


3.02 
4.35 
3.14 
3.40 
3.34 


0.031 
0.028 
0.242 
0.123 
0.070 


Fig. 3 
Fig. 3 



Table 1: Dependence of reconstruction performance on source plane pixel size, simulated ring 
image noise, psf width, and degree of regularization. 
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3.2. 



Fixed mass 



The main purpose of the simulations is to illustrate the linear inversion step, the 'inner 
cycle', i.e. that part of the semi-linear inversion method that differs from previous methods. 
For this reason in this sub-section we fix the mass parameters at the input values. The 
image inversion, therefore, is achieved in a single step using equations (5) and (12), for the 
unregularized and regularized cases respectively, and using equations (9) and (17) for the 
source covariance matrix. We consider the full problem, solving also for the mass parameters 



The unregularized inversion of the reference problem is provided in the remaining panels 
of Figure 1. The lower left panel shows the reconstructed source and the upper right panel 
shows the image of the reconstructed source convolved with the psf, i.e. the mm.—Xirn model 
fit to the simulated image. The bottom right panel shows the source a image i.e. the standard 
deviation in each pixel. This provides a visual impression of the uncertainties - note how 
the region of lowest a is bounded by the inner caustic. However the whole covariance matrix 
is required for a proper interpretation of the results. The lower middle panel is the source 
S/N image. In all the Figures 1 — 3, for the source a and S/N images the grayscale covers 
the full range of numbers in the panel. For the other panels the same grayscale range is used 
in each figure, to allow comparison of the relative noise levels. 

We measure several quantities to assess the quality of the inversion, listed in the re- 
maining columns of Table 1. The reduced x"^ in the image plane, xlmi^) provided in col. 
(7). The quoted uncertainty is given by \pljv where v is the number of degrees of freedom 
i.e. the number of image pixels (3626) minus the number of source pixels (225 or 900). 
The reduced in the source plane, Xso(^)> ^^d its uncertainty, is provided in col. (8). To 
account for the covariance terms this is computed using 



where Asj is the residual in the ith pixel. Here the number of degrees of freedom is the 
number of source pixels. Col. (9) provides the SjN summed over the small box in the 
source plane shown in Figure 1. The noise is computed as the square root of the sum of 
the elements in the covariance matrix, formed by stripping out from the source covariance 
matrix C the rows and columns corresponding to the pixels in the box. Col. (10) provides 
the absolute value of the significance As/cr of the worst-fit source pixel, and col. (11) lists 



in §3.3. 



3.2.1. Unregularized inversion 




(18) 
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the r.m.s. of the residuals in the source plane. 

The results for the reference problem, line (1) in Table 1, are all satisfactory: The 
reduced values in the image and source planes are both consistent with 1.0, and the 
significance of the worst pixel 2.80(7 (col. 10) is not unexpected given that there are 225 
source pixels. The summed S/N in the source box is an improvement on S/Nim- This might 
be expected since the box is restricted to the small region of the smirce plane containing 
nearly all the signal. At the same time it shows that the S/N is not greatly degraded by 
amplification of noise in the deconvolution step. We return to this issue below. We interpret 
these results as meaning that the inversion has succeeded and produced the correct solution 
to the well-posed problem of finding the source-pixel counts that give the best fit to the 
image. 

In simulation 2 we doubled the noise in the image plane. Comparing lines (1) and (2) 
in the table, the effect of this is to double the noise in the source plane (col. 11), and so 
halve the S/N of the detected source (col. 9), as expected. 
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X (arcsec) X (orcsec) X (arcsec) 

Fig. 2. — The plot shows the unregularized solution for the same problem as in Figure 1, but 
with source pixels half as large, and corresponds to line (5), Table 1. The source plane, top 
left panel and bottom row, is 0.75" x 0.75" with 0.025" pixels, and is centered on the optic 
axis. The source comprises two circular Gaussian components and is shown top left. Also 
marked is the line of the inner caustic for the isothermal ellipsoid lens. The image, convolved 
with the psf, FWHM 0.08", and with noise added, is shown upper middle. The image pixel 
size is 0.05" and the image box size is 5.0" x 5.0". The lower left panel is the source light 
distribution reconstructed from the image by semi-linear inversion without rcgularization. 
The grayscale range is the same as in Figure 1. The reconstruction is poor, because the 
source pixel size is too small. The upper right panel is the image of this source, convolved 
with the psf. The lower right panel displays the la uncertainty for the source pixels, and 
the lower middle panel is the source S/N image. The dotted square is the region over which 
S/Nso is measured. In each of Figures 1-3, counts in pixels in both the image and source 
plane are in units of surface brightness. 
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Variation of psf FWHM. We have investigated the effect of varying the width of the 
psf. Line (3) provides the results for no psf, and hne (4) provides the results for a psf 
FWHM of 0.16", double the reference value. Comparison of lines (1), (3), and (4) shows 
that as the psf FWHM increases the noise in the source plane, col. (11), increases and the 
source detection S/N, col. (9), decreases. This is as expected: In Fourier space the effect 
of the psf is to suppress the amplitude of the power spectrum of the source for large wave 
numbers. Therefore in the deconvolution process the noise on these scales is amplified. As 
the psf is broadened the power suppression is greater, and so the noise amplification in the 
deconvolution step is greater. The reduction in S/Ngo from 111.2 (hne 3) to 79.9 (hne 1), in 
going from no psf to psf FWHM of 0.08", is quite modest. This demonstrates that satisfactory 
inversion of ACS images using 0.05" source pixels is possible without regularization. 

Regardless of the degree of amplification of noise the various statistical quantities in 
cols (7)-(ll) of Table 1, lines (l)-(4), are all reasonable. This shows that in these cases 
the inversion is well behaved, and in none of the cases is the matrix F singular. This 
contrasts with the usual inversion problem, for example image deconvolution. With image 
deconvolution the number of parameters to solve for (the counts in the deconvolved image 
pixels) is typically the same as or greater than the number of constraints (the number of 
image pixels). In lensing, because of magnification, the number of image pixels may be much 
greater than the number of source pixels. This suggests, further, that in regions where the 
magnification is greatest it would be possible to use source pixels smaller than the image 
pixels. We consider this issue below. 

The results of these first four simulations indicate that, in some circumstances, provided 
the psf is not too broad, unregularized semi-linear inversion provides a useful solution. 

Variation of source pixel size. Figure 2 shows the same problem as Figure 1 but with 

0. 025" source pixels rather than 0.05" pixels. The results are summarized in line (5) of Table 

1. The quality of the reconstruction, lower left panel, is now dramatically worse, and outside 
the central region is clearly unsatisfactory. (The grayscale range of this panel is the same 
as in Figure 1.) Compared to line (1) the noise in the source (col. 11) has risen by a factor 
8, whereas intuitively one would expect only a factor 2 increase (4 times as many pixels). 
This is indicative of large amplification of noise, because the psf has sTipprcsscd the signal 
on these scales. This is a consequence of the fact that the separation in the image plane of 
the images of two adjacent source pixels is smaller than the psf size. Put another way, a 
resolution element in the image plane, traced back to the source plane, is oversampled by the 
source pixel size. The source covariance matrix now contains large, predominantly negative, 
covariance terms which correspond to the odd/even appearance in the outer regions of the 
source plane. 
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The high noise level in the outer parts of the source-plane belies the usefulness of this 
image. In fact the source is strongly detected, albeit at reduced S/Ngo = 20, even though not 
readily apparent to the eye. The source is clearly visible in the S/N image, however. At the 
same time the values in both the image and source planes remain satisfactory. Because 
of the larger magnification, the reconstruction is much better within the caustic line. This 
suggests it would be advantageous to use a variable pixel size across the source plane. For 
example, with reference to Figure 2, a scheme where the pixel size is 0.05" outside the caustic 
and 0.025" inside might be appropriate. We need to identify a criterion for choosing the pixel 
size that avoids the excessive amphfication of noise evident in Figure 2. There are clearly 
three variables which determine the minimum source pixel size: The image pixel size, a, the 
psf FWHM, b, and the magnification, c. We have had some success with a scheme which 
relates the source pixel size to the variable max{a,h/2)/c^l'^ . The results will be reported 
elsewhere (Dye and Warren, in prep.). 
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Fig. 3. — The plot shows regularized solutions for the same problem as in Figure 2, with different degrees 
of regularization. The middle row is for N\ = 1 (corresponding to line (8) of Table 1), and the bottom is for 
N\ = 3 (line (9) of Table 1). The source plane, top left panel, and middle and bottom rows, is 0.75" x 0.75" 
with 0.025" pixels, and is centered on the optic axis. The source comprises two circular Gaussian components 
and is shown top left. Also marked is the line of the inner caustic for the isothermal ellipsoid lens. The 
image, convolved with the psf, FWHM 0.08", and with noise added, is shown upper middle. The image pixel 
size is 0.05" and the image box size is 5.0" x 5.0". The left middle panel is the source light distribution 
reconstructed from the image by semi-linear inversion with regularization, A'^;^ = 1. The solution is much less 
noisy than the unregularized solution. Figure 2. The upper right panel is the image of this source, convolved 
with the psf. The middle right panel displays the la uncertainty for the source pixels. The center panel of 
the middle row is the source S/N image. The bottom row is the set of corresponding source-plane images 
for the case N\ — 3. Note the larger errors in the outermost band in the bottom right panel. This is a 
consequence of the choice of a gradient regularization term, since these pixels have fewer neighbours. The 
dotted square is the region over which S/Ngo is measured. In each of Figures 1-3, counts in pixels in both 
the image and source plane are in units of surface brightness. 
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3.2.2. Regularized inversion 

We now include linear regularization in the inversion. All the results reported here 
used the gradient form, equation (14). The results are quite similar for the different linear 
regularizing schemes described in §2.1.1, however. 

Lines (6) and (7) in Table 1 are the results for the reference problem with different 
degrees of regularization. As the regularizing term increases, the source becomes smoother 
and Xim increases. Line (6) is for Nx = I. Comparing line (6) to line (1) we see that the 
effect of regularization is to suppress the noise in the reconstructed source and to increase 
substantially the source S/N, col. (9). This is at the expense of a poorer match to the true 
source hght profile, as measured by cols (8) and (10). For — 2 the agreement with the 
input source is no longer acceptable. 

In hnes (8) to (10) we provide solutions for source pixel size 0.025", and psf FWHM of 
0.08", and different degrees of regularization, — 1, 3, 5. These results compare directly to 
the unregularized solution to the same problem, line (5). The solutions for simulation no. 
8, A^A = 1, and no. 9, A^^ = 3, are shown in Figure 3. The visual improvement, comparing 
the sequence of Figure 2 (unregularized). Figure 3 middle row (regularzsed, A^^ — 1), and 
Figure 3 bottom row (regularized, A^^ = 3), is dramatic. 

Comparing hnes (8) to (10) against line (5) we see that, again, regularization successfully 
suppresses noise, increasing the S/N of the detection of the source. As A^a increases, in this 
case xL increases only very slowly, much more slowly than in the case for larger pixels. This 
is partly due to the fact that we chose a smooth source, and the results would be different 
for a source with more small-scale structure. Nevertheless, it indicates that the standard 
criterion for the degree of regularization to apply, A^a = 1, is somewhat arbitrary. 

To summarise this sub-section, using a realistic problem, we have validated the theory 
of the linear inversion step set out in §2. This is the step that differs from the maximum- 
entropy method of Wallington et al. (1996), and therefore is the main point of the paper. 

3.3. Mass cycle 

In the present sub-section we report the results of solving the complete problem i.e. 
determining both the mass profile and the source light distribution. 

We first consider the unregularized case. Referring back to the example of Figure 1, the 
problem is to invert the image at upper middle. The free parameters are the five parameters 
describing the mass: x, y, ellipticity, position angle, and velocity dispersion. We searched 
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through the parameter space to find the min.— fit. At the minimum the matrix F supphes 
most of the terms of the curvature matrix. Following the precepts of §2.1.2, the remaining 
terms were filled in by measuring the relevant second partial derivatives of the surface. 
We found the surface to be completely smooth and parabolic near the minimum. The full 
curvature matrix was inverted to obtain the full covariance matrix for all the parameters, the 
mass terms as well as the counts in the source pixels. We found the input mass parameters 
were correctly recovered to within the uncertainties. We checked the full covariance matrix 
against the results of Monte-Carlo simulations and found excellent agreement. This confirms 
that, provided the chosen source pixel size is not too small, the unregularized semi-linear 
inversion method is a practical solution to the problem of inversion of a gravitational lens 
image with a resolved source. 

We also compared the terms in the source covariance matrix C = against the corre- 
sponding terms in the full covariance matrix. The differences are relatively small. Therefore 
the matrix C, at the global minimum, provides an approximation to the true source-pixel 
errors that may be very useful in the exploration stage, when considering different mass 
models and different pixelizations. 

In the regularized case we found, generally, that the procedure converged more rapidly 
than in the unregularized case. Regularized inversion can produce solutions which arc not 
true representations of the source (§3.4, Table 1). Nevertheless, we found, in contrast, that 
the solution for the mass parameters is very insensitive to the degree of rcgularization. In 
the regularized case, the curvature of the merit function cannot be used to obtain the full 
covariance matrix so that an alternative approach such as a Monte Carlo method must be 
adopted (§2.2.2). So the source covariance matrix, equation (17), at the global minimum, is 
particularly useful here as an approximation to the true source-pixel errors. 

3.4. Regularized vs unreguleirized 

We now include a debate on the relative advantages of the unregularized and regularized 
approaches. This may seem surprising given the excellent results achieved with rcgulariza- 
tion (comparing Figures 2 and 3). The weakness of the unregularized inversion is that in 
deconvolving the psf, the noise at large wavenumbers is boosted. The rcgularization term 
in the merit function imposes smoothness on the solution. In effect, the deconvolution (di- 
vision in Fourier space) is limited to the smaller wave numbers. However, this means that 
any real structure in the source at large wavenumbers is also suppressed. We are imposing 
a prejudice that the source is smooth and this might not be justified (see comments in §1). 
Rcgularization introduces -|-ve covariance between adjacent pixels, forcing the counts to be 
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similar. The regularized solution, then, is not so different to the unregularized solution with 
larger pixels. In this respect, it is interesting to compare the unregularized solution with 
0.05" source pixels (Figure 1), with the Nx = 1 regularized solution with 0.025" source pix- 
els (Figure 3). Noting that the two values of the source S/N are quite similar (79.9, 64.5, 
respectively. Table 1), it is a debatable point whether there is more information in the latter 
figure. 

A further point to note is that the regularized inversion can produce solutions which 
are satisfactory in terms of the fit to the image but which are not true representations of the 
source in the sense that the xio unsatisfactory. For example in line (7), Table 1, both xio 
and I As/a I are unsatisfactory. Measured by the same statistics, none of the unregularized 
inversions in the Table is unsatisfactory. The unregularized inversion gives a noisier but 
unbiased solution for the source light distribution, while the regularized inversion gives a 
smoother but biased solution. 

Despite regularization biasing the source, in the full mass cycle we find that the mini- 
mized mass parameters show little sensitivity to the degree of regularization. Furthermore, 
the regularized solution has the advantage that it converges more quickly. We discuss the 
practical significance of these two points in §4. Associated with this is the fact that regular- 
ization allows source pixel sizes of almost any size, unlike the unregularized case when pixel 
size must be chosen carefully. This can yield further speed advantages in the initial stages 
of an analysis, before the solution is refined. 

Overall we consider there are important advantages to using both regularized and un- 
regularized inversion in exploring the solution to a particular problem, and the choice will 

depend on the question being posed and any a priori knowledge concerning the source. Per- 
haps equally importantly, however, it makes sense to match the source pixel size to the data 
information content in terms of the S/N at different wavenumbers, or, in other words, to 
vary the pixel size depending on the magnification, as suggested in §3.2.1. 

4. Summary and recommendations 

We have developed a new method for the inversion of gravitationally-lensed images 
of extended sources for the case where the source light profile is pixelized. The method 
separates the linear dimensions of the problem (the counts in the source pixels) from the 
non-linear dimensions (the mass parameters). The method has been extended in a natural 
way to allow linear regularization of the inversion. The core of the routine is the procedure 
for inverting an image given a fixed mass profile. We have shown that this step, including 
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deconvolution of the psf, with or without regularization, is a hnear one. Since this step is 
usually achieved by searching the source parameter space for the merit-function minimum, 
the solution is reached much more quickly. The non-linear part of the problem has been 
reduced to the search for a minimum in the space of the mass parameters only. In the case of 
unregularized inversion, the full covariance matrix for all the (source+mass) parameters can 
be obtained very quickly. In the case of regularized inversion, a useful approximation to the 
covariance matrix for the source counts is obtained very simply, but Monte Carlo methods 
are needed to obtain the full covariance matrix. 

How the scmi-lincar method should be applied in practice depends on the problem 
posed. If one is interested in the quantitative details of the source light profile, for example, 
whether some apparent feature is real, then we recommend the unregularized solution. This 
is because regularization produces source profiles which are too smooth. Without regular- 
ization, an optimal source pixel size should be chosen. Having too large a source pixel may 
cause interesting detail to be lost. However, if the source pixel size is too small, the inverted 
image may have low S/N because of amplification of noise in the deconvolution step. If, 
on the other hand, one is interested only in the mass parameters, a regularized solution 
would be the appropriate choice: The mass parameters are rather insensitive to the degree 
of regularization and one benefits from an increase in inversion speed. 

Another consideration is that pixclizing the source uses a large number of parameters. 
As a rule, one is interested in finding the model with the smallest number of parameters 
that provides a satisfactory fit to an image. Therefore, in many cases, one might simply use 
the semi-linear method of inversion to provide an image of the source to guide the choice of 
parameterization. Here, again, the regularized solution might be the preferred option. 

In general, because it is so much easier to implement (§2.2.1), we recommend using 
zeroth-order regularization. Nevertheless, other considerations may override simplicity. The 
zeroth-order regularization term, in common with the maximum-entropy regularization 
term, is a local measure, independent of the counts in adjacent pixels. This can be an 
advantage or a disadvantage, depending on the actual light profile in the source. 

In the simulations presented here, we have used square source pixels which form a regu- 
lar grid. However, since the resolution across the image plane is fixed while the magnification 
varies, the resolution across the source varies. Therefore, to maximize the information con- 
tent in the reconstruction of the source it is necessary to use a variable source pixel size. We 
will present an analysis of semi-linear inversion with variable source pixel size in a future 
paper (Dye and Warren, in prep.). 

We have benefited from discussions with Paul Hewett, Geraint Lewis, Leon Lucy, and 
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Randall Wayth. 
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